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We study statistical properties of vibrated granular monolayers using molecular dynamics simu- 
lations. We show that at high excitation strengths, the system is in a gas state, particle motion is 
isotropic, and the velocity distributions are Gaussian. As the vibration strength is lowered the sys- 
tem's dimensionality is reduced from three to two. Below a critical excitation strength, a gas-cluster 
phase occurs, and the velocity distribution becomes bimodal. In this phase, the system consists of 
clusters of immobile particles arranged in close-packed hexagonal arrays, and gas particles whose 
energy equals the first excited state of an isolated particle on a vibrated plate. 

PACS: 81.05.Rm, 05.70.Fh, 02.70.Ns 



Granular media, i.e, ensembles of hard macroscopic 
particles exhibit rich, interesting, and only partiaUy un- 
derstood collective behavior ||^-^. The dynamics of 
driven or excited granular media is particularly impor- 
tant since the external energy source balances the energy 
loss due to collisions. Collective behavior of such systems 
is therefore important for establishing a more complete 
theoretical description of granular media. Here, we fo- 
cus on monolayer geometries where collapse, clusteriiig, 
and long range order have been observed recently [0-p|. 
In particular, experimental studies reported that the ve- 
locity distribution function may exhibit both Gaussian 
and non-Gaussian behavior under different driving con- 
ditions. 

In this study, we carry out molecular dynamics sim- 
ulations of vertically vibrated granular monolayers. To 
validate the simulation method, we verified the experi- 
mentally observed transition from a gas-like phase in high 
vibration strengths, to a cluster-gas phase in low vibra- 
tion strengths 0,^,0 • Additionally, we checked that sev- 
eral other details including the transition point, and the 
statistics of the horizontal velocities agree quantitatively 
with the experimental results. In particular, the hori- 
zontal energy vanishes linearly near the transition point, 
and the corresponding velocity distribution changes from 
a Gaussian to a non-Gaussian as the vibration strength 
is reduced. 

In contrast with the experimental studies, the simu- 
lations enables us to probe the vertical motion, an im- 
portant characteristic of the dynamics. Our results show 
that as the system approaches the transition point, the 
vertical energy drops by several orders of magnitude. 
Furthermore, the gas-cluster phase is characterized by 
a coexistence of clusters of immobile particles, and ener- 
getic gas particles, whose energy can be understood by 
considering an isolated particle on a vibrated plate. We 
also find that the deviation from the Gaussian behavior 
in the gas phase is directly related to the development 
of an anisotropy in the motion, i.e, significant differences 
between the horizontal and the vertical velocities. 

To study the dynamics of vibrated monolayers, we used 
the standard molecular dynamics simulation technique 



prt. We considered an ensemble of N identical weakly 
deformable spheres of mass m, radius R, and moment of 
inertia / ^ ^mB? . The simulation integrates the equa- 
tions of motion for the linear and angular momentums, 

™r« = T.j^t^'ij + "^gz, and /w, = T^j^-J^j x ^\j^_ re- 
spectively. Here, r^ is the position of the ith particle, 
LOi is its angular velocity and g is the gravitational ac- 
celeration. The force due to contact with the jth par- 
ticle in the direction normal (tangential) to the vector 
Vij = Yj — Vi is denoted by F^- (F^). The force be- 
tween two particles is nonzero only when they overlap, 
i.e., Fij = when |ry| > 2R. When there is an over- 
lap, the normal contact force F^- = FJ^*** -I- F^^"*^ be- 
tween the particles is a sum of the following forces: (a) 
A restoring force, FJ^''* = Ymi{\rij\ — 2i?)rjj/|r.y|, with 
Y the Young's modulus, and (b) An inelastic dissipa- 
diss _ „, ^.,,n 'Yhe tangential force is 



five force, Ylf"" ^ 
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the frictional force F*- = pshcar ^ — ^^miV*-. In the 
above, vf^- = {-Vij ■ ry)ry/|ry p, and v*^. = v^ - vf^ 
are projections of the relative velocities in the normal 
and tangential directions, respectively. The coefficients 
7„ and 7s account for the dissipation due to the relative 
motion in the normal and tangential directions, respec- 
tively. Overall, the molecular dynamics method has the 
advantage that it is amenable for parallel implementa- 
tion, and that it allows handling of collisions involving 
an arbitrary number of particles. 

Initially, particles are randomly distributed on the vi- 
brating plate, with a filling fraction p. The velocities were 
drawn independently from an isotropic Gaussian distri- 
bution. The plate undergoes harmonic oscillations in the 
vertical direction according to Zp{t) — A{t)cos{ujt) with 
A the vibration amplitude, and uj ~ I-kv with v the fre- 
quency. When a particle collides with the plate, it expe- 
riences the same force as if it were to collide with another 
particle moving with the plate velocity. The simulation 
was carried in a finite box with a height chosen to be 
large enough so that no collisions can occur with the box 
ceiling. Periodic boundary conditions were implemented 
horizontally. Overall, the simulation parameters were 
chosen to be as compatible as possible to the experimen- 



tal values [§: i? = 0.595mTO, g = g.Sm/s^, Y = lO^/s^, 
f = 70Hz, N = 2000, p = 0.463, 7^ = 100/s, and 
7„ = 200/s. The above normal dissipation parameter 
leads to a restitution coefficient oir — 0.95. Without loss 
of generality, we set the particle mass to unity, m = 1. 
We verified that the results reported in this paper were 
independent of the value of most of these parameters, as 
well as the nature of boundary conditions. 

We are primarily interested in statistical properties of 
the system in the steady state, and especially their de- 
pendence on the vibration strength, which can be quan- 
tified by the dimcnsionless acceleration F = Auj'^/g. The 
quantity F can be tuned by varying either to or A. We 
chose to fix the frequency and vary the vibration ampli- 
tude. This method should be valid as long as the time 
scale underlying the variation is larger than the systems' 
intrinsic relaxation time scales. Each of our simulations 
was initially run at ~ 10'^ oscillation cycles at a constant 
amplitude Aq. Then the amplitude was slowly reduced 
in a linear fashion according to A/Aq = 1 — t/r, with 
the decay time r = lO'^s (or alternatively r^ 10^ cycles). 
Throughout this paper we report measurements of aver- 
age quantities such as the temperature and the velocity 
distribution. These were obtained by averaging over 100 
consecutive oscillation cycles. 





Fig. 1. The gas versus the cluster-gas phase. Shown are 
the instantaneous particle positions (left) and the cumulative 
positions over 100 consecutive oscillation cycles (right) for vi- 
bration intensities F — 1.0 (up) and F = 0.6 (down). 

Qualitatively, we observe that above a critical vibra- 



tion intensity, F > Fc, particles are in a gas phase, in 
which their motion is random, as seen in Fig. 1. When 
F < Fc in addition to particles in the gas phase, hexago- 
nal ordered clusters form, as shown in a snapshot of the 
system. Furthermore, particles inside these hexagonal 
clusters are stationary, while particles outside the clus- 
ters move appreciably. This behavior is rather robust 
as it is independent of many of the underlying parame- 
ters including the dissipation parameters. Nevertheless, 
the simulations indicate that the critical acceleration Fc 
is primarily determined by 7„, while the stability of the 
clusters is governed by 7=,. Additionally, the relaxation 
time scale r had to be sufficiently small for the system 
to be able to fully relax any transient behavior. 
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Fig. 2 The vertical and horizontal temperatures versus the 
vibration strength. 

Experimental measurements of the horizontal temper- 



ature, defined by Th 



indicate a linear de- 



pendence in the vicinity of the critical point M 
rH(x(F-Fc). 



(1) 



Our simulations confirm this linear behavior, as shown 
in Fig. 2 and Fig. 3. Furthermore, the horizontal energy 
practically vanishes below the transition point, as this 
quantity decreases by 3 orders of magnitude for F < Fc. 
This is reminiscent of a sharp phase transition and it is 
therefore sensible to view Fc as a critical point. This lin- 
ear behavior can be used to estimate the critical point, 
and a linear least-square-fit yields Fc = 0.763, in good 
agreement with the experimental observation, Fc — 0.77 
[Q . We conclude that the near critical behavior observed 
numerically agrees both qualitatively and quantitatively 
with the experimental observations. 

Simulations also allow measurements of the vertical ve- 
locities. We find that the vertical energy Ty = '2{vz'^) 
decreases sharply near the transition point as well. How- 
ever, in contrast with the horizontal energy, it does not 
vanish below the transition point. Therefore, the veloc- 
ities develop a strong anisotropy as the vertical and the 



horizontal velocities behave quite differently. This re- 
flects the fact that the system is far from equilibrium. 
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Fig. 3. Near critical behavior of the horizontal tem- 
perature. The critical acceleration was determined to be 
Tc = 0.763 from the linear fit Th = B{T - T^), plotted as 
a solid line. 

More detailed velocity statistics is provided by the ve- 
locity distribution. We observe that at high accelerations 
the particle motion is nearly isotropic, i.e., the ratio of 
horizontal to vertical energy is of the order unity. In- 
deed, this ratio approaches a value of roughly 0.5 (see 
Fig. 2). In addition, the velocity distribution is Gaus- 
sian (see Fig. 4) and the system is practically three- 
dimensional. However, as the acceleration is decreased, 
the two-dimensional geometry becomes more and more 
pronounced. The vertical motion dominates over the hor- 
izontal one, and the horizontal velocity distribution de- 
parts strongly from a Gaussian distribution. Near the 
phase transition point the large velocity tail becomes 
nearly exponential (see Fig. 4). Below the transition 
point, a significant fraction of the particles have a nearly 
vanishing horizontal velocity, and the distribution of hor- 
izontal velocities is strongly enhanced near the origin. 

The deviation from the Gaussian behavior can be 
quantified using the kurtosis, defined via the fourth and 
second moments of the distribution, k = (u^)/(w^)^. In- 
deed, in the limit of high vibration intensities, F ^ 1, 
this parameter approaches the Gaussian value k — )■ 3. 
On the other hand, near the phase transition point, 
i.e, as F — > Fc, this parameter approaches the expo- 
nential value K — > 6. It proves useful to examine how 
the kurtosis depends on Th/Tv, the ratio between the 
horizontal and the vertical energies. As shown in the 
inset to Fig. 4, the smaller the ratio (or equivalently, 
the larger the anisotropy), the larger the deviation from 
a Gaussian distribution. Hence, whether the velocity 
distribution is Gaussian or not reflects the degree of 
anisotropy in the particle motion. Non-Gaussian distri- 



butions has been observed experimentall 
retically [|ll]-|3|, and numerically §,P-[l5| 
three-dimensional geometries. 
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Fig. 4 The distribution of horizontal velocities in the gas 
state, and in the vicinity of the critical point. The velocities 
were normalized by the RMS velocity Wrma = {v^)^ ■ A Gaus- 
sian distribution is plotted as a reference. The inset plots k 
the kurtosis of the distribution of horizontal velocities versus 
Th/Tv the ratio of horizontal to vertical energies. The two 
extremal points correspond to the data sets plotted in this 
figure, i.e, V = 0.8 and F = 6. 
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Fig. 5 The distribution of vertical velocities in the cluster- 
gas phase. 

The distribution of vertical velocities can be used to 
distinguish between cluster and gas particles. Indeed, the 
vertical velocity distribution changes its character from 
a unimodal to a bimodal distribution in the cluster-gas 



phase (see Fig. 5) with the low velocity peak correspond- 
ing to the cluster particles, and the high velocity peak 
corresponding to the gas particles. Interestingly, the lo- 
cation of the high velocity peak does not change as F 
decreases. In fact, this energy can be understood by con- 
sidering the first excited state energy of a single particle 
bouncing on a vibrating plate, which can be calculated 
to be § 



-■4(^r 



(2) 



or in our case Ei = 8.16 cm? / s^ . We verified this result 
by simulating the motion of a single particle on a vibrat- 
ing plate. Interestingly, the energy of the gas particles 
falls within less then 10% of this value. Therefore, below 
the phase transition point particles residing in clusters 
are in the ground state, i.e., they are moving with the 
plate. Furthermore, the rest of the particles constituting 
the gas phase are in the first excited state of an isolated 
ball on a vibrating surface. This indicates that particles 
in the gas phase are essentially noninteracting. 

The vertical velocity distribution can be used to study 
the fraction of particles in each phase by simply integrat- 
ing the area under the respective energy peaks. As shown 
in Fig. 6, Pq, the fraction of particles in the gas phase, is 
almost independent of the vibration intensity below the 
transition point. As the transition point is approached 
this fraction rapidly decreases and ultimately vanishes 
for F » Fc. Although this quantity does not undergo a 
sharp transition, its behavior is consistent with our pre- 
vious estimate of the transition point from the horizontal 
energy behavior. 
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Fig. 6 The fraction of particles moving with the plate. The 
fraction Po was obtained from the horizontal velocity distri- 
bution by integrating the area enclosed under the low velocity 
peak (see Fig. 5). 

In summary, we have studied the dynamics of vibrated 
granular monolayers using molecular dynamics simula- 



tions. We find that the transition between the gas and 
the cluster-gas phases can be regarded as a sharp phase 
transition, and that the horizontal energy decreases lin- 
early near the transition point. We have shown that at 
high vibration strengths, the particle motion is isotropic, 
and the velocity distributions are Gaussian. The devia- 
tion from a Gaussian distribution were found to be closely 
related to the degree of anisotropy in the motion. We 
have also shown that below the phase transition, the ve- 
locity distribution is bimodal. The cluster particles move 
with the plate, while the gas particles behave as a nonin- 
teracting gas, as their energy agrees with the first excited 
state of an isolated vibrated particle. 

Our results agree both qualitatively and quantitatively 
with the experimental data. This shows that the underly- 
ing phenomena can be explained solely by the simulated 
interactions, i.e., contact force interactions, no attrac- 
tive forces, and dissipative collisions. Other mechanisms, 
possibly present in the experiment, such as electrostatic 
forces, etc. are therefore not responsible for the phase 
transition. It will be interesting to use molecular dy- 
namics simulations to determine the full phase diagram 
of this system by varying the density and the driving 
frequency. 
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